arXiv: 1508.00817vl [physics.med-ph] 4 Aug 2015 


1 


Joint Reconstruction of Absorbed Optical Energy 
Density and Sound Speed Distributions in 
Photoacoustic Computed Tomography: 

A Numerical Investigation 

Chao Huang, Kun Wang, Member, IEEE, Robert W. Schoonover, Lihong V. Wang, Eellow, IEEE, 

and Mark A. Anastasio, Senior Member, IEEE 


Abstract —Photoacoustic computed tomography (PACT) is a 
rapidiy emerging bioimaging modaiity that seeks to reconstruct 
an estimate of the absorbed opticai energy density within 
an object. Conventionai PACT image reconstruction methods 
assume a constant speed-of-sound (SOS), which can resuit in 
image artifacts when acoustic aberratious are significant. It has 
been demonstrated that incorporating knowiedge of an object’s 
SOS distribution into a PACT image reconstruction method 
can improve image quaiity. However, in many cases, the SOS 
distribution cannot be accurateiy and/or convenientiy estimated 
prior to the PACT experiment. Because variations in the SOS 
distribution induce aberrations in the measured photoacoustic 
wavefields, certain information regarding an object’s SOS dis¬ 
tribution is encoded in the PACT measurement data. Based on 
this observation, a joint reconstruction (JR) problem has been 
proposed in which the SOS distribution is concurrently estimated 
along with the sought-after absorbed optical energy density from 
the photoacoustic measurement data. A broad understanding 
of the extent to which the JR problem can be accurately and 
reliably solved has not been reported. In this work, a series of 
numerical experiments is described that elucidate some important 
properties of the JR problem that pertain to its practical 
feasibility. To accomplish this, an optimization-based formulation 
of the JR problem is developed that yields a non-linear iterative 
algorithm that altematlngly updates the two image estimates. 
Heuristic analytic insights into the reconstruction problem are 
also provided. These results confirm the ill-conditioned nature 
of the joint reconstruction problem that will present significant 
challenges for practical applications. 

Index Terms —Photoacoustic computed tomography, optoa- 
coustic tomography, ultrasound tomography, image reconstruc¬ 
tion 


I. Introduction 

Photoacoustic computed tomography (PACT), also known 
as optoacoustic or thermoacoustic tomography, is a rapidly 
emerging hybrid imaging modality that combines optical im¬ 
age contrast with ultrasound detection Q-lH. In PACT, the 
to-be-imaged object is illuminated with a short laser pulse that 
results in the generation of internal acoustic wavefields via the 
photoacoustic effect. The amplitudes of the induced acoustic 
wavefields are proportional to the spatially variant absorbed 
optical energy density distribution within the object, which 
will be denoted by the object function A(r). The acoustic 
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waves propagate out of the object and are subsequently mea¬ 
sured by use of wide-hand ultrasonic transducers. The goal 
of image reconstruction in PACT is to estimate the object 
function A(r) from these measurements. 

Image reconstruction methods for PACT are often based 
on idealized imaging models that assume an acoustically ho¬ 
mogeneous medium ii-ilHi. However, these assumptions are 
not warranted in certain biomedical applications of PACT 0. 
Numerous image reconstruction methods have been proposed 
mni-iini that compensate for aberrations of the measured 
photoacoustic (PA) wavefields caused by an object’s speed- 
of-sound (SOS) variations, c(r), and hence improve PACT 
image quality. It has been demonstrated that these methods 
can improve the fidelity of reconstructed images by incorpo¬ 
rating accurate knowledge of the SOS variations in the PACT 
imaging model. However, accurate estimation of c(r) prior to 
the PACT study generally requires solution of an ultrasound 
computed tomography (USCT) inverse problem, which can 
present experimental and computational challenges lfT8l - ll2^ . 

An important observation is that, because variations in the 
SOS distribution induce the PA wavefield aberrations, certain 
information regarding an object’s SOS distribution is encoded 
in the PACT measurement data. Based on this observation, 
it is natural to question whether A{r) and c(r) can both 
be accurately determined from only the PACT measurement 
data. Il2^ - ll27l . thereby circumventing the need to perform a 
dedicated USCT study. This will be referred to as the joint 
reconstruction (JR) problem and is the subject of this article. 

Theoretical and computational studies of the JR problem 
have been conducted but all are limited in scope. Theoretical 
work on the JR problem that neglects discrete sampling effects 
has established that A{r) and c(r) can be uniquely determined 
from the measured PACT data only under certain restrictive 
assumptions regarding the forms of A(r) and c(r) lfT4l . Il28l 
or the measurement surface ll29l . However, the uniqueness of 
the JR problem for the general case has not been established. 
Another study established that the solution of the linearized 
JR problem is unstable ll^ and suggested that the same 
conclusion would hold for the general case where wavefield 
propagation modeling is based on the full wave equation. 

Despite the lack of theoretical works, others have moved 
forward and developed computational methods for solving 
the JR problem by use of discretely sampled measurement 
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data Il2^ - ll25l . In ||2^ . an iterative reconstruction method was 
proposed to jointly estimate both A(r) and c(r). That study 
employed a geometrical acoustics propagation model and 
assumed a priori information regarding the singular support 
of c(r). In ll24l . ||25l, a JR method based on the Helmholtz 
equation was proposed that was solved by the finite element 
method (FEM). While this method is grounded in an accurate 
model of the imaging physics, it suffers from an intensive 
computational burden. A similar JR approach was proposed 
E6i that employed a time-reversal (TR) adjoint method. All 
of these works are preliminary in the sense that they did 
not systematically explore the numerical properties of the JR 
problem or provide broad insights that allow one to predict 
when accurate JR may be possible. In combination with the 
scarcity of theoretical works, this indicates an important need 
to further elucidate the practical feasibility of JR. 

To address this, the primary objective of this work is to 
investigate the numerical properties of the JR problem, which 
will provide important insights into its practical feasibility. A 
novel JR method is developed for this purpose. The developed 
reconstruction method is based on an alternating optimization 
scheme, where A(r) is reconstructed by use of a previously- 
developed full-wave iterative method while c(r) is recon¬ 
structed by use of a nonlinear optimization algorithm based on 
the Frechet derivative of an objective function with respect to 
c(r) ED, E2- Computer-simulation studies are conducted to 
investigate the topology of the cost function defined in the 
optimization-based approach to JR. Additionally, numerical 
experiments are conducted that reveal how the supports and 
relative smoothness of A(r) and c(r) affect the numerical 
stability of the JR problem. Demonstrations of how errors 
in the imaging model associated with imperfect transducer 
modeling and acoustic attenuation affect JR accuracy are also 
provided. 

The paper is organized as follows. In Section [III the 
imaging physics of PACT in acoustically heterogeneous media 
is reviewed briefly. The derivation of the Frechet derivative 
with respect to c(r) of a pertinent objective function is also 
provided. Section |III] describes the alternating optimization 
approach for solving the JR problem. In Sec. HYl heuristic 
insights into how the relative extents of the spatial supports 
of A(r) and c(r) can affect the ability to perform accurate 
JR are provided. The computer-simulation methodology and 
numerical studies are given in Secs. IV] and |V|] The paper 
concludes with a summary and discussion in Section IVTIl 


II. Background 


A. Photoacoustic wavefield propagation in heterogeneous me¬ 
dia 


We consider PA wavefield propagation in lossless fluid 
media having a constant mass density. Fet p{r, t) denote 
the photoacoustically-induced pressure wavefield at location 
r G and time t > 0. The photoacoustic wavefield p{r, t) 
satisfies ID: 


V2p(r,f) 


1 d'^p{r,t) 

c(r)2 


( 1 ) 


subject to initial conditions 

p(r,0) = r(r)A(r), 


dp{r, t) 


dt 


= 0 , 


( 2 ) 


t=o 


where r(r) is the Grueneisen parameter that is assumed to be 
known. 


B. Frechet derivative with respect to c(r) 

Here, for simplicity, we neglect the acousto-electrical im¬ 
pulse response (FIR) and the spatial impulse response (SIR) of 
the ultrasonic transducers employed to record the PA signals. 
However, the impact of these will be addressed in Section 
IVI-Dl The quantity p(r™,f) represents the PA data recorded 
by the m-th transducer at location r™ (m = 1, • • ■ ,M). For 
ease of description, we represent the measured PA data as con¬ 
tinuous functions of t, but the results below will be discretized 
for numerical implementation as described in Section |III] 

The problem of reconstructing c(r) from PA data for a fixed 
A{r), defined as Sub-Problem #2 below, can be formulated 
as an optimization problem in which the following objective 
functional is minimized with respect to c(r); 


M 

^[c(r)] = / 

m=l "'0 




( 3 ) 


subject to the constraint that p{r'^,t) satisfies Eqs. ID and 
@1, where T denotes the maximum time at which the PA data 
were recorded. 

Gradient-based algorithms can be utilized to minimize the 
nonlinear functional E)- Such methods require the functional 
gradient, or Frechet derivative, of £ with respect to c(r), which 
can be calculated by use of an adjoint method ED, OH- In 
the adjoint method, the adjoint wave equation is defined as 




1 9^g(r, f) 

dt^ 




subject to terminal conditions 


g(r,T) = 0, 


dq{r,t) 

dt 


= 0 . 

t=T 


( 4 ) 

( 5 ) 


The source term s{r,t) is defined as 


M 


s(r, f) = > ^ [p(r-, f) - p(r™, t)]S{r - r™). 


m—1 


(6) 


Upon solving ID and 01, the Frechet derivative of £ with 
respect to c(r) can be determined as lED, ED, 


Vr.£ = - 


c(r)3 


^^dp{r,t) dq{r,t) 


( 7 ) 


dt dt 

Once the Frechet derivative is obtained, it can be utilized by 
any gradient-based method as the search direction to iteratively 
reduce the functional value of 0 . The computation of this 
Frechet derivative from discrete measurements is described in 
the Appendix. 


III. Joint Reconstruction of A(r) and c(r) 

In this section, a JR method for concurrently estimating 
A(r) and c(r) is formulated based on an alternating optimiza¬ 
tion strategy. 
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A. Discrete imaging model 

Let the A X 1 vectors 

A = [A(ri),-- 

• ,Al(r^r 

(8) 

and 

c= [c(ri),-- 


(9) 


denote the finite-dimensional approximations of A(r) and c(r) 
formed by sampling the functions at the N vertices on a 
Cartesian grid that correspond to locations r„, n = 1, • • • ,N. 
As introducted earlier, r^, fc = 1, • • • , M denote the transducer 
locations. 

The quantity 

Pi = (10) 

represents the measured PA data sampled at time t = lAt 
(Z = 1, • • • , L) at each transducer location. Here, At is the 
sampling time step, and L is the total number of time steps. 
The complete set of measured PA data can be represented by 
the LM X 1 vector 

p= [pi,--- (11) 

By use of ([8]l, (|9]l, and (fTTT) . a discrete PACT imaging model 
can be expressed as na 

p==H(c)A, (12) 

where H(c) is the LM x N system matrix that depends non- 

linearly on c. A procedure to establish an explicit matrix 
representation of H(c) was provided in Ifia . 

In conventional applications of PACT the SOS distribution 
c is assumed to be known. Alternatively, the goal of the 
JR problem is to concurrently estimate A and c from the 
measured data p by use of the model in Eq. (fT2l) . 

B. Optimization-based joint image reconstruction 

Based on Eq. (fTSI) . the JR problem can be formulated as 

A,c = argmin ||H(c)A-p|P + AiRa(A) + A 2 Rc(c), (13) 

A>0,c>0 

where Ra(A) and Rc(c) are penalty functions that impose 
regularity on the estimates of A and c, respectively, and 
Ai, A 2 are the corresponding regularization parameters. As 
discussed in Sec. IVI-Al the cost function in ( fT3l l is non- 
convex. However, a heuristic alternating optimization approach 
can be employed to find solutions that approximately satisfy 
(foi l. This approach consists of the two sub-problems described 
below: (1) reconstruction of A given c, and (2) reconstruction 
of c given A. 

Sub-Problem #1: Reconstruction of A given c: The 
problem of estimating A for a given (i.e., fixed) c can be 
formulated as the penalized least squares problem 

A = argmin ||H(c)A - p|p -f AaRa(A), (14) 
A>0 

where Aa is the regularization parameter, which is different 
from Ai in Eq. ( fOl ) in general. Reconstruction methods have 
been proposed for solving problems of this form m, El. 


Sub-Problem #2: Reconstruction of c given A: Eor a 

given A, an estimate of c can be formed as 

c = argmin ||H(c)A - p||^ -|- AcRc(c), (15) 
00 

where Ac is the regularization parameter, which is different 
from A 2 in Eq. (fOl) . in general. Equation (fTSl l can be solved by 
use of gradient-based methods, which require computation of 
the gradient of the objective function in Eq. (ffSl l with respect 
to c. Details regarding this gradient computation are provided 
in Sec. III-BI and the Appendix. 


Algorithm 1 Alternating optimization approach to JR of A 
and c_ 

Input: p, A(°\ ca, Cc, Aa, Ac 

Output: A, c 

1: i = 0 

2 : while ei < ca and £2 < Ec do 

3: A(*+^^ •(—Fa(A(A, {Sub-Problem #1} 

4: ^ Fc(c('\ a('+^\ p, Ac) {Sub-Problem #2} 

5: ei ^Dist(A«,A(‘+i)) 

6: £2 ^Dist(c«,c(‘+1)) 

7: i ^ I -I- 1 

8 : end while 
9: A ^ A^®) 

10: C •(— 


Alternating optimization algorithm: As described in Algo¬ 
rithm 1, JR of A and c can be accomplished by alternately 
solving Sub-Problems #1 and #2. The quantities A*^°^ and 
are the initial estimates of A and c, respectively, and 
£A and £c are convergence tolerances. The functions ‘Fa’ and 
‘Fc’ compute the solutions of Eqs. (fT4l) and ( fTSl ), respectively, 
and are described below in Section [V] The function ‘Dist’ 
measures the Euclidean distance between A*^®) and A*^®+^^ (or 
between c^®^ and c^®+^)). 

IV. Heuristic insights into support conjecture 
AND Sub-problem #2 

As revealed throughout this work, difficulities in solving 
Sub-Problem #2 represents a significant challenge for JR. In 
this section, heuristic insights into how the relative extents of 
the spatial supports of A(r) and c(r) can affect the ability to 
accurately solve Sub-Problem #2, and hence the JR problem, 
are provided. 

The supports of the functions A(r) and c(r) will be denoted 
as supp(A) and supp(c), respectively. The supports are defined 
to be the regions where A{t) ^ 0 and c(r) — cq 0, 
respectively. Here cq is the known SOS in the background 
(i.e., in the coupling water-bath). The functions A(r) and 
c(r) are assumed to be compactly supported, reflecting the 
physical constraint that their extents are limited and the object 
of interest resides within the imaging system. 

The following assumptions will be made in order to permit 
a simple analysis of the problem that yield insights. Eirst, it is 
assumed that acoustic scattering is weak. More specifically, a 
geometrical acoustics model Eol is employed to describe the 











4 


propagation of the PA wavefields. Amplitude variations are 
assumed to be negligible and wavefront aberrations are mod¬ 
eled by computing the time-of-flight along straight propagation 
paths. The second assumption exploits the fact that an arbitrary 
function A(r) can be decomposed into a collection of point 
sources. Specifically, it is assumed that the PA signal generated 
by each of these point sources can be recorded independently 
by transducers. In reality, of course, this is not the case. This 
assumption, in effect, assumes that an ‘oracle’ records the PA 
signals and is able to decompose them into the individual 
components that were produced by each point source. Below, 
we will describe how the relative extents of supp(A) and 
supp(c) affect the ability of the oracle to accurately solve Sub- 
Problem #2. Since the oracle has access to more information 
than is actually recorded in an experiment, inability of the 
oracle to perform accurate image reconstruction implies that 
accurate image reconstruction by use of the actual recorded 
data will also be unfeasible. This is the logical basis for the 
analysis below. Third, the analysis is presented in 2D but can 
be extended to the 3D case readily. Without loss of generality, 
the measurement surface is assumed to be a circle with radius 
R that encloses supp(A). 

For convenience, we introduce the perturbed slowness dis¬ 
tribution s(r) = - -As assumed above, the oracle can 

resolve the PA signal generated by each point source that 
comprises A{r). As such, the time it takes for a PA pulse 
to propagate from its emission location G supp(A) to each 
transducer location rp = [i?cos(/3), i?sin(/3)]'^ (/3 G [—7r,7r)) 
is also known to the oracle, where the geometry is defined in 
Fig.m From these time-of-flight (TOF) data, a tomographic 
data function can be defined as 

3(^/3,= Tcoiri3,rf3,) - T{rs,rf3) - T{rs,r'p), (16) 

where is the transducer location defined by the intersection 
of the line connecting rp and r^, and the measurement circle. 
The quantity Tca{vp^Ypi) denotes the time it takes for a 
pulse to propagate between the two transducer locations when 
only the background medium is present. Since a straight-ray, 
geometrical acoustics, model is assumed, this data function is 
related to the sought-after slowness distribution as 

3(i‘/3,r^;i’s) = / s(r)dr, (17) 

where the path of integration is along the line segment 
L{r'p,rp) that connects the two transducer locations. 

Two cases will be considered that reveal the impact of the 
relative extents of supp(A) and supp(c) on the ability of the 
oracle to reconstruct s(r), or equivalently, c(r), given A(r). 
The first case corresponds to the situation where supp(c) C 
supp(A), as depicted in Fig. [T] In this case, for example, 
supp(A) could correspond to the area occupied by a soft 
tissue structure in which the SOS is approximately the same 
as the background SOS and supp(c) would correspond to a 
region within the tissue possessing a different SOS. It can 
be verified readily that for every r G supp(c), a collection 
of Tg G supp(A) exists such that values of the data function 
that specify the line integrals along all paths L(r'p,rp) that 
intersect r are accessible to the oracle. Stated otherwise, in 



Fig. 1. Schematic of the 2D circular measurement geometry employed in the 
heuristic analysis of Sub-Problem #2 described in Sec. hy] The coordinates 
rj 3 and denote the transducer locations that correspond to the intersection 
points of the line r • h = d with the measurement circle. In this example, 
supp(c) C supp(A). 

this case, all projection data that are needed to uniquely invert 
the 2D Radon transform m in Eq. dnli are accessible and 
therefore s(r) or, equivalently, c(r) can be exactly determined 
in a mathematical sense. 

In fact, the requirement supp(c) C supp(A) can be relaxed 
to supp(c) being enclosed by supp(A), as shown in Fig. |2ta). 
This is because we only require G c)supp(A), where 
i9supp(A) denotes the boundary of supp(A), to determine the 
complete set of projection data as described above. 


(a) (b) 

Fig. 2. Two cases addressed in the heuristic analysis of Sub-Problem #2 
described in Sec. m (a) Case where supp(c) is enclosed by supp(A); and 
(b) Case where supp(c) is not enclosed by supp(A). In (b), certain line integral 
data for the slowness distribution in a subset of supp(c) are not measured, as 
indicated by the dashed lines through the regions covered by lines. 

The second case corresponds to the situation where supp(c) 
is not enclosed by supp(A), as shown in Fig. |2|b). In this 
case, it can be verified that the complete data function cannot 
be determined, indicating that the values of all line integrals 
through s(r) are not accessible as depicted in Fig. |2b). That 
example is analogous to the interior problem in X-ray CT, 
which possesses no unique solution ll34l . 

These observations are summarized as follows. 

Support conjecture:: When geometrical acoustics is valid, 
Sub-Problem #2 cannot be accurately solved when supp(c) is 
not enclosed by supp(A). 

V. Computer-Simulation Studies 

A general description of the computer-simulation method¬ 
ologies is described below. Although the optimization ap¬ 
proach to JR described above is based on the 3D wave 
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equation, 3D JR remains computationally intensive. For com¬ 
putational convenience, the 2D formulation is investigated in 
this work. The specific numerical experiments designed to 
investigate the numerical properties of the JR problem are 
described subsequently in Sec. [VT] 


Xeon(R) E5645 CPUs and an NVIDIA Tesla C2075 graphics 
processing unit (GPU). The GPU was equipped with 448 
1.15 GHz CUDA cores and 5 GB global memory. The Jacket 
toolbox was employed to accelerate the computation of 
([T]i and (|4]i on the GPU. 


A. Simulation of idealized PA measurement data 

Numerical phantoms were utilized to represent A(r) and 
c(r), which were described by Eqs. ([8]) and (|9]l. These phan¬ 
toms, denoted by A and c, contained 512 x 512 pixels with 
a pitch of 0.25 mm. (i.e., N = 512^ in Eqs. ([8]l and ®.) 
Eor a given A and c, the lossless PA wave equation was 
solved by use of the MATLAB k-Wave toolbox iTSl to produce 
simulated PA signals at each transducer location. Each PA 
signal contained 6000 temporal samples recorded with a time 
step Af = 50 ns. The measurement geometry consisted of 
800 transducers that were evenly distributed on the perimeter 
of a square with a side length of 100 mm. The object 
was contained within this region. Note that, unless stated 
otherwise, the generation of the simulated PA measurement 
data in this way avoided an ‘inverse crime’, due to a different 
choice of discretization parameters from those employed by 
the reconstruction algorithm as described below. 

B. Simulation of non-idealized PA measurement data 

To investigate the properties of JR under more realistic con¬ 
ditions, additional PA data sets were simulated that considered 
certain physical factors. Measurement noise was modeled by 
adding 3% (with respect to the maximum value of the noiseless 
data) white Gaussian noise (AWGN) to the simulated PA data. 
Additional factors considered included acoustic attenuation, 
the effect of the electrical and spatial impulse responses of 
the transducers. Specific details are provided in Sec. IVI-DI 


VI. Numerical investigations and analyses 
A. Topology of the JR cost function in Eq. ( 1731 ) 

The effectiveness of the optimization-based approach to JR 
depends on the topology of the cost function that is minimized 
in Eq. (fl^ . In these studies, only the data fidelity term was 
considered (Ai = X 2 = 0). If the cost function is not convex or 
quasi-convex, a global minimum (i.e., an accurate JR solution) 
may not be returned by the optimization algorithm when the 
initial estimate is not close enough to the global minimizer. 
Moreover, when the cost function is not strictly convex, there 
is no guarantee that the solution is unique. In this section, 
a low-dimensional stylized example is considered to yield 
insights into the general characteristics of the cost function 
topology. 

The numerical phantoms shown in Eig. [3 were employed 
to represent A and c. To establish a low-dimensional rep¬ 
resentation of these phantoms, a discretization scheme that 
differed from Eqs. ® and (|9]l was employed (in this study 
only). Namely, A was described as A = where Ag 

represented the value within the uniform disk and Af, is the 
known (assumed to be zero) background value. Similarly, c 
was described as c = [cg, Cb] where Cg represented the value of 
the uniform annulus and Cb is the fixed background value out¬ 
side of the annulus. Eor each (Ag,Cs) pair, simulated PA data 
p were computed. Subsequently, the value of ||H(c)A — p|p 
was computed and plotted as a function of the scale factors 
of Ag and Cg. 


C. Implementation details for image reconstruction 

The implementations of the two image reconstruction sub¬ 
problems in Algorithm 1 are described below. The function 
‘Fc’ that computes the solution of (fTSI l was implemented based 
on the MATLAB k-Wave toolbox llT5l . Specifically, the wave 
equation ([T]i and the adjoint wave equation (|4]l were solved 
numerically by use of the k-space pseudospectral method. 
The computed PA wavefield and the adjoint wavefield were 
employed to compute the gradient of the objective function in 
(fTSl l (see Appendix). The gradient was subsequently utilized 
by the limited-memory BEGS (L-BEGS) algorithm to solve 
([I5]i Ebl-iEsl. The implementation of the function ‘Fa’ that 
solves (fTrli can be found in ifT^ . In this study, a total variation 
(TV) penalty was adopted. The ‘Dist’ function measured the 
difference in terms of root mean squared error (RMSE), and 
the convergence tolerances ca and Cc were empirically chosen 
to have a value of 10“^ throughout the studies. In all studies, 
the initial estimates of A and c were set to be = 0 and 
c(o) = 1480 m/s, which is the background SOS. Both A and 
c were reconstructed on a uniform grid of 256 x 256 pixels 
with a pitch of 0.5 mm. 

All simulations were computed in the MATLAB environ¬ 
ment on a workstation that contained dual hexa-core Intel(R) 


o 


'1500 

(a) (b) 

Fig. 3. Subfigures (a) and (b) display the numerical phantoms that specified 
A and c that were employed to investigate the JR cost function topology in 
Sec. rVLAl 

Eigure HJa) displays a surface plot of the data fidelity term 
||H(c)A — p|p as a function of the scale factors of Ag and Cg, 
while profiles corresponding to Ag = 1 are shown in subfigure 
(b). The results show that the data fidelity term is not convex 
with respect to Cg. These results suggest the cost function may 
also be non-convex with respect to more general and higher¬ 
dimensional SOS distributions c. 

To investigate how the spatial structure of A can influence 
the topology of the cost function, the above procedure was 
repeated when the radius of the disk in A was varied. Profiles 
of the normalized data fidelity term corresponding to Ag = 1 
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Fig. 4. Subfigure (a) displays a surface plot of the data fidelity term in the 
JR cost function and the corresponding profile corresponding to As = 1 is 
shown in subfigure (b). Details are provided in Sec. IVI-AI 


are displayed in Fig. Ha). These data reveal that the ‘valley’ 
containing the global minimum of the cost function widens as 
the radius of the disk heterogeneity in A increases. This can be 
seen more clearly in the plots of the width of the valley versus 
the disk radius that are shown in Fig. Eh). This observation 
suggests that the quality of the initial estimate for c can be 
relaxed when the radius of the support of A(r) increases. 




Fig. 5. Profiles of the normalized data fidelity term corresponding to As = 
1 for different radii of the disk phantom that specified A are displayed in 
subfigure (a). Subfigure (b) displays a plot of the width of the valleys in (a) 
versus the disk radius. 


In summary, these results establish that, for a fixed A, 
the JR cost function is generally non-convex with respect to 
c. Accordingly, Sub-Problem #2 is generally a non-convex 
problem. Sub-Problem #1 is convex if the penalty is convex. 
The overall non-convexity of the JR problem indicates that 
accurate initial estimates of c will generally be necessary 
to avoid local minima that represent inaccurate solutions. 
Additionally, the non-convexity of the problem suggests that 
uniqueness of the solution is not guaranteed. 

B. Numerical investigations of Sub-Problem #2 

As described above, Sub-Problem #2 corresponds to a non- 
convex optimization problem that can be difficult to solve 
in practice. Accordingly, errors that arise when solving Sub- 
Problem #2 can accumulate in Algorithm 1 and hinder the 
ability to perform accurate JR. Below, numerical experiments 
are reported that reveal insights into some mathematical prop¬ 
erties of A(r) and c(r) that influence the ability to accurately 
solve Sub-Problem #2. 

1) Effect of spatial supports of A(r) and c(r).' Studies 
were conducted to investigate the extent to which the support 
conjecture, provided in Sec. HYl influences the ability to ac¬ 
curately solve Sub-Problem #2, and hence perform JR, by use 
of perfect measurements. The numerical phantom employed 


to represent the SOS distribution, c, is shown in Fig. |6] The 
SOS values were representative of human breast tissues. We 
first considered two choices for A, shown in Figs. |7ja) and 
|71d), which were designed to satisfy the support conjecture 
(supp(c) is enclosed by supp(A)). For each choice of A, a 
corresponding estimate of c was estimated by solving Sub- 
Problem #2 from noiseless simulated PACT measurements. 
The value of the regularization parameter Ac in Eq. (fTSI l was 
zero. The reconstructed estimates of c corresponding to the 
two choices of A are shown in Figs. |7jb) and [Tjc). The 
corresponding profiles extracted from the central rows of the 
reconstructed estimates are displayed in Figs.|7jc) and|7|f), re¬ 
spectively. These results reveal that it is possible to reconstruct 
accurate estimates of c from perfect PACT measurements via 
Sub-Problem #2 when the specified A satisfies the support 
conjecture. The effects of noise and other measurement errors 
are addressed later. 



Fig. 6. The numerical SOS phantom described in Sec. IVI-Bl that was utilized 
in the investigations of Sub-Problem #2. 
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Fig. 7. Numerical investigations of Sub-Problem #2 - Spatial support effects: 
Two different phantoms A that satisfy the support conjecture ai‘e shown in 
subfigures (a) and (d). The estimates of c obtained by solving Sub-Problem 
#2 by use of noiseless data and the corresponding A are shown in subfigures 

(b) and (e). Image profiles thi'ough the estimates of c are shown in subfigures 

(c) and (f). 


We next considered the four choices for A, shown in the left 
column of Fig. [8] (Figs. [Sja), (d), (g), and (j)), which, roughly 
speaking, were designed to violate the support conjecture 
to different extents. As above, corresponding unregularized 
estimates of c were estimated by solving Sub-Problem #2 from 
noiseless simulated PACT measurements. The reconstructed 
estimates of c corresponding to the different choices of A 
are shown in the middle column of Fig. The corresponding 
profiles extracted from the central rows of the reconstructed 
estimates are displayed in the right column. Figures |3^b) 
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and (e) reveal that, despite that the violation of the support 
conjecture, accurate estimates of c could still be reconstructed. 
Figures [8jh) and (k) demonstrate that the accuracy of the 
reconstructed estimate of c degrades as the size of the support 
of A(r) is reduced. 
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Fig. 8. Numerical investigations of Sub-Problem #2 - Spatial support effects: 
Four different phantoms A that violated the support conjecture are shown in 
subfigures (a), (d), (g), and (j). The estimates of c obtained by solving Sub- 
Problem #2 by use of noiseless data and the corresponding A are shown in 
subfigures (b), (e), (h), and (k). Image profiles through the estimates of c are 
shown in subfigures (c), (f), (i), and (1). 


In summary, these observations reveal that, in general, the 
support conjecture does not need to be satisfied in order to 
achieve accurate reconstruction of c for an exactly known 
A and perfect measurement data. This may be explained 
by the fact that the support conjecture was established by 
use of geometrical acoustics that represents a simplification 
of the complicated acoustic wave propagation in both our 
simulations or in practice. However, in the cases where the 
support conjecture was satisfied, c was accurately estimated. 
Moreover, the simulation results indicate that the extent to 
which the supports of A(r) and c(r) overlap does affect the 
ability to accurately solve Sub-Problem #2, and hence the JR 
problem. 

2) Effect of relative spatial bandwidths of A{r) and c{r): 
Studies were conducted to investigate the extent to which the 
relative spatial bandwidths of A(r) and c(r) influence the 
ability to accurately solve Sub-Problem #2 by use of perfect 
measurements. Figure |9] shows the numerical phantoms for A 
and c. To exclude the effects related to the supports of A(r) 
and c(r), the phantoms were designed to satisfy the support 
conjecture. The spatial structures of the original phantoms 


in Fig. |9] are identical, indicating that their spatial frequency 
bandwidths are identical. The phantom depicting A was subse¬ 
quently convolved with different Gaussian kernels to generate 
additional phantoms that possessed different spatial frequency 
bandwidths. The full width at half maximum (FWHM) of the 
blurring kernel was employed as a summary measure of the 
relative spatial bandwidths of the smoothed phantoms. Sub- 
Problem #2 was subsequently solved with Ac = 0 by use of 
both perfect and noisy simulated measurement data for each 
of the smoothed versions of A. 

2400 
2200 
2000 
1800 
1600 

(a) (b) 

Fig. 9. The numerical phantoms representing (a) A and (b) c that are 
described in Sec. IVI-B2I 



The reconstructed estimates of c are shown in Fig.[T0l From 
top to bottom, the results correspond to relative bandwidth 
ratios (of A to c) of 0.25, 0.44, and 1.0, respectively. Figures 
Ma), (e), and (i) correspond to images reconstructed from 
perfect measurement data and the associated image profiles are 
displayed in Figs. lTOt b). (f), and (j). Figures fTOl' cl. (g), and (k) 
correspond to images reconstructed from noisy measurement 
data and the associated image profiles are displayed in Figs. 
[Tol' dl. (h), and (1). The RMSE of the reconstructed c with 
respect to the relative bandwidth ratio of A to c for both 
noiseless and noisy cases is displayed in Fig. [TT] 
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Fig. 11. Numerical investigations of Sub-Problem #2 - Spatial bandwidth 
effects: Plots of RMSE of the estimated c as a function of the bandwidth 
ratio of A to c. Subfigures (a) and (b) correspond to the noiseless and noisy 
results, repectively. 


These results reveal that the relative spatial frequency band- 
widths of A(r) and c(r) affect the numerical stability of Sub- 
Problem #2 and, hence, that of the JR problem. The accuracy 
of the reconstructed c in the present of measurement noise 
for a specified A was observed to be influenced strongly by 
the relative spatial bandwidths of A and c. Specifically, in 
order to accurately reconstruct c in the presence of noise, the 
presented results suggest that the spatial frequency bandwidth 
of A(r) should be comparable to or larger than the spatial 
frequency bandwidth of c(r). We will refer to this as the k- 
space conjecture hereafter. 
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Fig. 10. Numerical investigations of Sub-Problem #2 - Spatial bandwidth effects: As described in Sec. |VI-B2| unregularized estimates of c were reconstnacted 
from noiseless and noisy data and a series of A that had different spatial bandwidths relative to that of the sought-after c. From the top to the bottom rows, 
the results correspond to relative bandwidths of A to c of 0.25, 0.44, and 1.0, respectively. Subfigures (a), (e), and (i) correspond to images reconstructed 
from perfect measurement data and the associated image profiles are displayed in subfigures (b), (f), and (j). Subfigures (c), (g), and (k) con'espond to images 
reconstructed from noisy measurement data and the associated image profiles are displayed in subfigures (d), (h), and (1). 


3) Effect of perturbations of A{r): In the studies described 
above, perfect knowledge of A was assumed. Below, a numer¬ 
ical experiment is described that provides insights into how 
small perturbations in the assumed A affects the accuracy of 
the reconstructed c obtained by solving Sub-Problem #2. 


Figures WB a) and (c) display two similar numerical phan¬ 
toms depicting A. The RMSE between these phantoms is 
0.004. Perfect PACT measurements were simulated for each 
of the two A, for a given c (not shown). Figures fTSl b) and (d) 
display the reconstructed estimates of c when the A specified 
in Fig. fT2l' a) and (c) was assumed, respectively. These results 
demonstrate that the problem of reconstructing c for a given 
A is ill-conditioned in the sense that small changes in A can 
produce significant changes in the reconstructed estimate of 
c. This observation is consistent with the theoretical results in 
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It is interesting to note that the (A, c) pair in the top row 
of Fig. [reproduces nearly identical PA data, at all transducer 
locations, to that produced by the (A, c) pair in the bottom 
row. The simulated noiseless pressure data at an arbitrary 
transducer location produced by the two (A, c) pairs is shown 
in Fig. [13] where the pressure signals are observed to overlap 
almost completely. The RMSE between the two sets of PA 
data was RMSE = 3.2 x 10“"^. These results suggest that 
the solutions of the JR problem in PACT may not be unique. 
Consequently, it indicates that accurate JR of A and c, in 
general, may not be possible. 


Fig. 12. Numerical investigations of Sub-Problem #2 - Effect of perturbation 
of A: Two numerical phantoms representing A are shown in subfigures (a) 
and (c). As described in Sec. IVI-B3I these two phantoms are very similar, 
with a RMSE between them of only 0.004. Unregularized estimates of c 
reconstructed by use of noiseless simulated measurement data and the A 
specified in (a) and (c) are shown in subfigures (b) and (d), respectively. 
Although the phantoms depicting A are very similar, the reconstructed 
estimates of c are not. 


C. Feasibility of JR with idealized data 

The numerical instability of Sub-Problem #2, as examined 
in the previous section, implies that the solution to the JR 
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Fig. 14. Images obtained via JR from noiseless PA measurements: The first and third columns display the estimates of A and c, and the corresponding 
image profiles are displayed in the second and fourth columns, respectively. Each row corresponds to use of different regularization parameters and Ac. 
From top to bottom, (A^, Ac) are (0,10“®), (10“^, 10“®), (0,10“'*), and (10“®, 10“^), respectively. 



(a) (b) 


Fig. 13. Numerical evidence of non-uniqueness of the JR problem: Simulated 
PA measurement data were computed from the (A, c) pairs shown in Figs. 
Illj al and (b) and Figs. lllj cl and (d). The two pressure profiles con'esponding 
to an arbitrary transducer location are superimposed in subfigure (a). Subfigure 
(b) displays a zoomed-in version of subfigure (a). Similar agreement between 
the profiles was observed at all transducer locations. 

problem is also numerically unstable. This was confirmed by 
computing solutions to the JR problem from perfect measure¬ 
ment data. The data were perfect in the sense that they did 
not contain measurement noise. Moreover, ‘inverse crime’ was 
committed in which the same forward model and discretization 
parameters (pixel size =0.5 mm) were employed to produce 


the simulation data and to conduct image reconstruction. These 
data were produced by use of the phantoms shown in Fig. |9] 
which satisfied both the support and the k-space conjectures. 
Unregularized JR was performed by use of Algorithm 1 with 
Aa = Ac = 0 and the results are displayed in the top row 
of Fig. [14] Despite the use of perfect measurement data and 
a favorable choice of A and c, neither A nor c could be 
accurately reconstructed. 

Additional studies were conducted in which noiseless mea¬ 
surement data were computed without committing inverse 
crime. Namely, the pixel size employed to simulate the mea¬ 
surement data was 0.5 times the pixel size employed in the 
reconstruction algorithm. From these data, regularized esti¬ 
mates of A and c were computed. In rows 2-4 of Fig. [141 the 
corresponding regularization parameters are Aa = Ac = 10“®, 
10“^^, and 10“^, respectively. These results demonstrate that 
the numerical instability of the JR problem can be mitigated by 
incorporating appropriate regularization. However, the salient 
observation here is that regularization was required to obtain 
accurate image estimates, even if no stochastic measurement 
noise or other significant modeling errors other than discretiza- 
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tion effects were introduced. 

D. Feasibility of JR with imperfect data 

1) Effect of stochastic measurement noise: JR was per¬ 
formed by use of noisy versions of the PA data. The first 
and second rows of Fig. [15] display the estimates of A 
and c, respectively, along with corresponding image profiles. 
These results were obtained with regularization parameters 
A ,4 = 10“^ and Ac = 10“^. There results show that, with 
appropriate regularization, JR can be robust to the effects of 
AWGN. 



Fig. 15. Images obtained via JR from noisy PA measurements: The top and 
middle rows display the reconstructed estimates of A and c, respectively, 
along with the corresponding image profiles. The bottom row displays an 
estimate of A reconstnacted by a PACT reconstruction method that employed 
a constant SOS that was manually tuned to minimize RMSE. 


To compare with the JR results, A was also reconstructed 
by use of a full-wave PACT reconstruction method Uhl that 
neglected acoustic heterogeneity by employing a constant SOS 
value. This estimate of A is displayed in the third row of 
Fig [15] The reconstruction method assumed a constant SOS 
value of 1600 m/s, which was manually tuned to minimize 
the RMSE of the reconstructed image. The RMSE of the 
reconstructed A by use of the JR method and the PACT 
method ignoring SOS variations was 0.01 and 0.21, respec¬ 
tively. These results confirm that the estimate of A obtained 
by performing JR is more accurate than the corresponding 
estimate reconstructed by use of a PACT method assuming a 
constant SOS. 

2) Effect of other modeling errors: In practice, besides 
stochastic measurement noise due to the electronics, other 
forms of measurement data inconsistency will be present. Eor 


example, the imaging model utilized in this work neglects 
the spatial impulse response (SIR) of the transducers and 
acoustic attenuation within the object and coupling medium 
lEa, m-M- Additionally, because the pressure data are 
assumed to represent the measurable quantity, the electrical 
impulse response (EIR) of the transducer is assumed to be 
known exactly and the deconvolution process to obtain the 
pressure data is assumed to be error-free. These assumptions 
will be violated in a real-world experiment. A study was 
conducted to demonstrate the performance of JR based on 
our idealized imaging model in the presence of stochastic 
measurement noise and these factors that are neglected or only 
partially compensated for. 

The simulated PA data containing the effects of these 
physical factors was generated as follows. 

1) Simulated PA data were generated in a lossy medium. 
The A and c phantoms shown in Pig. [9]were employed. 
Acoustic attenuation was introduced by use of an acous¬ 
tic attenuation coefficient a that was described by a 
frequency power law of the form a(r,/) = 

ll45l . The frequency-independent attenuation coefficient 
at) = 10 dB MHz and the power law exponent y = 2.0 
were employed in the data generation, which correspond 
to the values of ao and y in human kidneys that have the 
strongest acoustic attenuation among typical biological 
tissues ll46l . The simulated PA data were contaminated 
by 3% AWGN. 

2) To model the effects of the SIR, the simulated PA 
data, described above, were computed on a grid with 
a pitch of 0.1 mm and recorded by 4000 transducers 
that were evenly distributed on the sides of a square 
with side length 100 mm. The recorded data from 
every 20 consecutive transducers were then averaged to 
emulate the SIR of a 2 mm line transducer. Note that 
this simplified model of the SIR accounts only for the 
averaging of the pressure data over the active surface 
of a transducer El and neglects other factors that may 
contribute to the SIR of a real-world transducer. 

3) The simulated PA data containing the acoustic atten¬ 
uation effects, the SIR effects and measurement noise 
were convolved with an EIR of an actual transducer 
iol, El- The convolved data were then deconvolved 
by use of a curvelet deconvolution technique B9l . In 
the deconvolution, a perturbed EIR was employed that 
was produced by adding 2% Gaussian noise into the 
spectrum of the original EIR. Pigure [161 a) displays the 
perturbed and original EIR and Pig. [Tbl bl displays the 
deconvolved and true pressure signal for a particular 
transducer location. 

Prom these data, JR of A and c was performed with 
regularization parameters A^ = 10“^ and Ac = 10“^. The 
JR results are displayed in the top and middle rows of Pig. 
[TT] These results suggest that, even with a sufficient A that 
satisfies the support and k-space conjectures, accurate JR may 
not be feasible in practice due to its instability unless model 
errors are small. However, the jointly reconstructed A has 
smaller RMSE = 0.12 compared to the iterative result (the 






































11 




Fig. 16. Panel (a): inaccurate EIR compared to the original EIR. Panel (b): 
deconvolved pressure data by use of the inaccurate EIR compared to original 
pressure data. 


bottom row of Fig. [TtT i that was reconstructed with constant 
SOS of 1600 m/s and has RSME = 0.22. This shows that, 
even though accurate JR may not be feasible in practice, the 
JR method provides the opportunity to improve the accuracy 
of the reconstructed A as compared to the use of a PACT 
reconstruction method that assumes a constant SOS EO). 



Fig. 17. Images obtained via JR from noisy PA measurements in the presence 
of model error: The top and middle rows display the reconstmcted estimates 
of A and c, respectively, along with the corresponding image profiles. The 
bottom row displays an estimate of A reconstructed by a PACT reconstruction 
method that employed a constant SOS that was manually tuned to minimize 
RMSE. 


VII. Conclusion and discussion 

Because variations in the SOS distribution induce aber¬ 
rations in the measured PA wavefields, certain information 
regarding an object’s SOS distribution is encoded in the 
PACT measurement data. As such, several investigators have 
proposed a JR problem in which the SOS distribution is 
concurrently estimated along with the sought-after absorbed 


optical energy density. The purpose of this work was to 
contribute to a broader understanding of the extent to which 
this problem can be accurately and reliably solved under 
realistic conditions. This was accomplished by conducting a 
series of numerical experiments that elucidated some impor¬ 
tant numerical properties of the JR problem. These studies 
demonstrated the numerical instability of the JR problem. 

The presented findings are consistent with and corroborate 
previous theoretical studies of the JR problem. Namely, Ste- 
fanov et al proved the instability of the linearized JR problem, 
which suggested instability of the more general JR problem as 
well lEa. In ED, a condition for unique reconstruction of c(r) 
given A(r) was provided, which is consistent with our support 
condition (see Theorem 3.3 therein). In previous work regard¬ 
ing the development of JR algorithms, Chen et al proposed a 
similar optimization-based approach to JR ||26l. They solved 
the optimization problem by use of an optimization algorithm 
called the TR adjoint method. Although their algorithm was 
different from ours, they obtained similar results; Accurate 
JR images were not produced when A is deficient, but the 
jointly reconstructed A could be more accurate than the one 
reconstructed by use of the TR method with a constant SOS. 

Similar and results can be found in the works by Jiang et 
al ll24ll . II 25 I . in which the authors proposed an optimization 
approach to JR that was based on the Helmholtz equation 
instead of the wave equation. By use of that method, the au¬ 
thors observed that the accuracy of JR results was affected by 
the temporal frequency band employed in the reconstruction. 
Specifically, the frequency ranges covering lower frequencies 
gave more accurate JR results than higher frequencies. This 
observation is implicitly contained in our heuristic k-space 
conjecture. Their observations and our k-space conjecture can 
be understood by noting that band-pass or high-pass filtered A 
are not physical because the non-negativity of A does not hold 
in those cases. Both results showed that the accuracy of JR is 
impacted by the spatial spectrum of A. When A and c possess 
sharp boundaries and the same structures (i.e., are simply 
scaled versions of each other), the authors also showed that 
biased estimates of A and c could be jointly reconstructed by 
incorporating Marquardt and Tikhonov regularizations into the 
reconstruction method. Note that in this case, A and c satisfied 
our support and k-space conjectures. By use of regularization, 
the authors showed that their algorithm was insensitive to ran¬ 
dom noise in the measurement, which is congruous with our 
observations. Although the reconstructed images were biased, 
the authors showed that the jointly reconstructed A was more 
accurate than the image reconstructed with a homogeneous 
SOS, which is consistent with our results. In addition, they also 
observed that the jointly reconstructed A was more accurate 
than the jointly reconstructed c, which, again, indicated the 
inverse problem of reconstructing c is more unstable compared 
to the reconstruction of A. 

However, none of these works reported a systematic in¬ 
vestigation of the numerical properties of the JR problem 
or provided broad insights that allow one to predict when 
accurate JR may be possible. In particular, the impact of the 
spatial support and spatial frequency content of A(r) relative 
to that of c(r) was not explored. In the current study, we have 
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demonstrated that, even if the measurement data are perfect, 
accurate JR may not be achievable. 

The investigation of the JR problem by use of experimental 
data remains a topic for future investigation. However, based 
on the presented studies, the task of performing accurate JR 
under experimental conditions is likely to be highly challeng¬ 
ing and will require accurate modeling of the imaging operator. 
This will generally require the use of the 3D wave equation 
instead of 2D wave equation. A line search is inevitable in 
any nonlinear optimization algorithm that is employed to solve 
Sub-Problem #2, which will create a very large computational 
burden in the 3D case. Additionally, in this study, a lossless 
fluid medium was assumed by the reconstruction method. 
However, in certain applications, density variations and/or 
acoustic absorption may not be negligible m, 1521. These 
assumptions can, in principle, be relaxed when formulating 
the JR problem. 

Finally, due to its instability, it will likely be beneficial 
to incorporate additional information into the JR problem. 
One possibility is to augment the PACT measurement data 
with a small number of ultrasound computed tomography 
(USCT) measurements. The investigation of the JR problem 
by combining PACT and USCT measurements is underway 

Ell. 
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Appendix; calculating the gradient of ( fTSl i 

The gradient of the first term in Eq. ( fTSl ) can be calculated 
by discretizing the the Erechet derivative (|7]) 

a||H(c)A-p|p ^ _4r;-3 o f V P^+i ~ P^-i 

dc \ 2 

/ = 1 

qi+i-qz-i , . s / X 

o-^-f (pi - Po) o (qi - qo) 

+(pl-i - pl-2) o( qi-i - qL-2)} 

where o denotes Hadamard product, is defined as 

C~^ = ,c{rN)~^r, (19) 

Pi and q; (/ = 0, • • • ,L — 1) are defined as 

Pi = [p{ri,lAt), ■ ■ ■ ,p(r7v, (20) 

and 

qi = [q{ri,lAt), ■ ■ • ,q{rN,lAt)]'^, (21) 

representing the PA wavefield and the adjoint wavefield sam¬ 
pled at the 3D Cartesian grid vertices r„ (n = 1, • • • , N) and 
at time t = lAt, respectively. 


If TV-penalty is adopted, the gradient of the second term in 
(m is given by ll5^ 


9Ac|c|tv 

dc 


— Ac [ci, 




( 22 ) 


and 

3 3 

c„ = (3[c]„ - ^[c]„-){e + ^([c]„ - [c]„-)2}-^ 

3 3 (23) 

-E + 

i=l j=l ^ 

where e is a small positive number to prevent the denominators 
being zeros, and [c]„ denotes the n-th grid node of c, and 
[c]^- and [c]^+ are neighboring nodes before and after the 
n-th node along the /-th dimension (i = 1,2,3), respectively. 
Likewise, [c], +,- denotes the neighboring node that is after 

)j 

the n-th node along the i-th dimension and before the n-th 
node along the j-th dimension. 

The gradient of the objective function in (fTSl l is then given 
by the sum of Eqs. (fTSl l and (l22li . 
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